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When classical systems fail to explore their entire configurational space, intriguing macroscopic phenomena 
like aging and glass formation may emerge. Also closed quanto-mechanical systems may stop wandering 
freely around the whole Hilbert space, even if they are initially prepared into a macroscopically large 
combination of eigenstates. Here, we report numerical evidences that the dynamics of strongly interacting 
lattice bosons driven sufficiently far from equilibrium can be trapped into extremely long-lived 
inhomogeneous metastable states. The slowing down of incoherent density excitations above a threshold 
energy, much reminiscent of a dynamical arrest on the verge of a glass transition, is identified as the key 
feature of this phenomenon. We argue that the resulting long-lived inhomogeneities are responsible for the 
lack of thermalization observed in large systems. Such a rich phenomenology could be experimentally 
uncovered upon probing the out- of- equilibrium dynamics of conveniently prepared quantum states of 
trapped cold atoms which we hereby suggest. 

The ergodicity axiom in classical statistical mechanics states that, during its time evolution, a closed mac- 
roscopic system uniformly explores the entire phase space compatible with conservation laws, so that the 
time average of any observable comes to coincide with the micro-canonical ensemble average and, when the 
observable is local, also with the canonical Gibbs ensemble average. Nonetheless, ergodicity can be violated in 
classical systems, a noticeable example being glasses 1 . Quantum effects might also spoil ergodicity by preventing 
the wave function from diffusing within all available configurations. This phenomenon is actually known to occur 
in the presence of disorder and manifests itself either by single-particle 2 or many-particle 3 " 5 wave function 
localization. However, alike classical models for glassy behavior, ergodicity breakdown in the quantum dynamics 
may not necessarily require disorder and it could instead be entirely due to frustrating dynamical constraints 6,7 . 
This issue is currently attracting great interest 8,9 , since well controlled realizations of closed quantum systems have 
become feasible upon trapping cold atomic species 10 . Indeed, similarly to what can be done in numerical 
simulations, one can prepare atoms in a given initial state and probe their time evolution under a 
Hamiltonian whose parameters are fully under control, thus offering the unique opportunity to monitor the 
ergodicity principle at work in the quantum realm 11 . 

In this work, we report numerical evidences that an isolated system of strongly interacting bosons, modeling 
atoms in optical lattices, can be trapped during its evolution into long-lived inhomogeneous metastable states, 
provided that its internal energy exceeds a certain threshold. We argue that the slowing down of high-energy 
incoherent excitations in the strongly correlated system is the key feature responsible for this dynamical arrest, 
much resembling a kind of glass transition. By formulating the problem in a different language, we explicitly show 
that a system initially prepared in a inhomogeneous state is unable to diffuse within the entire configurational 
space; such a dynamical localization in the many-body Hilbert space looks intriguing and may represent a kind of 
many-body Anderson localization 3 that occurs without disorder. The above phenomenon is put in further 
relation with and deemed responsible for the lack of ergodicity observed in large finite size systems. This belief 
is confirmed by means of a novel time- dependent variational Monte Carlo method that we introduce hereby and 
an experimental set up to uncover this very rich phenomenology is also suggested. 

One of the simplest models that can be realized in experiments is the Bose-Hubbard Hamiltonizan 10,12 : 

n= -/£ (^ + h.c.) + ^„,.(n,-l), (l) 

(Uj) i 

characterized by the amplitude /for a bosonic species of an atom to hop between nearest-neighboring wells of an 
optical lattice and by a local repulsion [/among atoms localized in the same potential well. The operators bj and b t 
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create and destroy, respectively, a boson on site i, and Yi[ = b\bi is the 
density operator 36 . Experiments are often performed with aniso- 
tropic lattices that realize a collection of almost uncoupled chains, 
a fortunate case for numerical simulations that we shall mainly con- 
sider hereafter, apart from a brief excursion in two dimensions 
towards the end of the paper. 

In one dimension, obstacles to ergodicity can arise in integrable 
models 13 . However, the Hamiltonian (1) is not integrable and, 
indeed, there are experimental evidences that its dynamical evolution 
may succeed in fast relaxing to a thermal state. Specifically, in a re- 
cent experiment a system of 87 Rb atoms, well described by the 
Hamiltonian (1), has been prepared in a state in which the sites of 
the optical lattice were alternatively empty and singly occupied 14 . 
This state was let evolve for different experimental conditions cor- 
responding to different ratios U/J. Even at the largest value U/J ~ 10, 
the initial density profile (... 1,0,1,0,. . .) was found to rapidly relax to 
the homogeneous thermal one (...|,^,^,|,...), with half a boson 
per site, much faster than the integrable counterparts of non-inter- 
acting or infinitely- interacting (i.e., hard-core) bosons and consis- 
tently with the increased number of relaxation channels that opens 
once integrability is lost. 

Numerical simulations of the above experiment successfully 
reproduce the observed thermal behavior 14 . However, there are also 
numerical evidences pointing out a breakdown of ergodicity in the 
same model (1) but within a different region of the parameter space, 
specifically when the number of bosons is one per site. This case, as 
well as any other one with integer density, is special because, at 
equilibrium and at zero temperature, the model (1) undergoes a 
quantum phase transition into a Mott insulator above a critical Ul 
J. Even though the phase transition is washed out by thermal fluctua- 
tions, nonetheless its influence on the spectrum seems to prevent 
thermalization above a certain LV/in numerical simulations of finite 
size systems 26,29 ; a result that may 27 or may not 28 prelude to a true 
breakdown of ergodicity in the actual thermodynamic limit. 

Here, we propose an experiment a lot alike the one previously 
described 14 , which may distinguish very sharply a change from cha- 
otic to non-chaotic dynamics in the Bose-Hubbard model (1). Our 
proposal starts by observing that, in the gapless phase next to the 
Mott insulator at one boson per site, low-energy itinerant Bogoliubov 
quasiparticles should coexist with high-energy incoherent excita- 
tions, which, for U / J^>1, can be identified as sites occupied by more 
than a single boson. It is well possible that the overall relaxation 
depends critically upon the population of those high-energy excita- 
tions 20 , which can be assessed by tailoring initial states with lots of 
doubly occupied sites rather than none as in the experiment of Ref. 14 . 
This is indeed what we propose and simulate by means of different 
and complementary numerical tools, such as exact diagonalization, 
time-evolving block decimation 15 and a novel time-dependent vari- 
ational Monte Carlo algorithm, as discussed in more detail in the 
Supplementary Material. 

Results 

Inhomogeneous initial states and dynamical localization. We start 
from analyzing the model at density n = 1, when at equilibrium a 
Mott transition occurs at a critical (U"//) c ~3.5 16 . We imagine to 
prepare an initial state where all the sites are either empty or 
doubly occupied. In particular we consider the states depicted in 
Fig. 1, namely with clusters of doubly occupied sites of variable 
size S d . These states are let evolve with the spatially homogeneous 
Hamiltonian (1) for different U/J, below and above the critical value. 
While for small U/J the density profile rapidly reaches the 
equilibrium configuration (...1,1,1,1,...), for large U/J, it stays close 
to its initial value for a remarkably long time. Eventually, since the 
system is finite, the density profile approaches the homogeneous 
plateau, with small residual oscillations that get damped as the 
system size increases. 



S d = l 



S d = 3 



Figure 1 | Inhomogeneous initial states constituted by clusters of doubly 
occupied sites (large blue circles) and empty sites (small dots). The size of 
each cluster is denoted by S d . 

We can define a relaxation time t r (whose inverse is shown in 
Fig. 2, for S d = 1), as the first time for which the local density 
approaches its homogeneous value. We highlight that z R , at a specific 
(U /J) d / n > has a sudden step up, which becomes sharper and sharper 
as the system size increases. The above results show that above 
(U /J) d / n the system has the tendency to stay dynamically trapped 
into long-lived inhomogeneous configurations. 

We emphasize that such a behavior is all the more remarkable not 
because doubly occupied sites seem unable to decay, which is known 
to require for U/J^>\ very long times and large system sizes 18 , but 
rather because they are not capable to move, hence restore trans- 
lational symmetry. We can better understand this surprising result 
by an effective Hamiltonian that can be derived following the same 
reasoning of Refs. 19 ' 20 . For a sufficiently large interaction, it is jus- 
tified to project the evolution onto states with the same potential 
energy per site (7/2, at least for time scales shorter than (T 2 // 3 . One 
realizes that states with the same potential energy but with triply and 
singly occupied sites start to contribute only at order /YU 3 , so that, 
with accuracy T 2 / (7, the number of doubly occupied sites is conserved. 
If we associate a fictitious spin up or down to a doublon (doubly 
occupied site) or a holon (empty site), respectively, we find that the 
effective Hamiltonian that controls the evolution reads: 



(ij) 



(2) 



which describes a hard-axis ferromagnetic Heisenberg model 19 . We 
note that the evolution of an XXZ spin chain starting from an anti- 
ferromagnetic ordered initial state has been studied numerically in 
Ref. 21 with DMRG. Here we consider the effective dynamics (2) in a 
regime when the anisotropy would favor a ferromagnetic ground 
state and when the initial state contains clusters of up and down 
spins of increasing lenght. The smaller cluster size of Sd = 1 corre- 
sponds to a Neel initial state considered in Ref 21 . In the right panel of 
Fig. 2 we show the time-dependence of the clusters density in the 
large U regime, evolved according to the effective Hamiltonian (2). 
Remarkably, we see that even for small clusters the system fails to 
restore the spatial homogeneity up to very long time scales, which 
turn to be far beyond those currently accessible in typical experi- 
mental setups. 

The slowing down of the dynamics can be traced back to the 
effective attraction among doublons 20 , which makes their aggregates 
hard to break up. In other words, what seems to matter more is the 
dissociation of clusters of doublons, rather than the decay of a single 
one. Indeed we have checked (not shown) that in the large U regime 
the system has the tendency to get dynamically stuck into clusters of 
doublons of finite size, whereas a much faster annihilation and 
recombination rate is observed in the small U limit. 

To get further insights into the dynamical behavior of the system, 
we recast the problem in a different language. Starting from the initial 
state, denoted as |0), we can generate an orthogonal basis set | /), i = 0, 
1,..., by repeatedly applying the Hamiltonian (see Supplementary 
Material). In this Lanczos basis of many-body wave functions, the 
Hamiltonian has the form of a tight-binding model on a semi-infinite 
chain. Each site i = 0, 1,... corresponds to a many-body state, it has 
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Figure 2 | Left panel - Inverse relaxation times t# 1 of the local density for the initial state (. . . 2, 0, 2, 0, . . .). Exact diagonalization results are reported, 
with darker points marking larger systems, respectively N = 8, 10 and 12 with periodic-boundary conditions. In the Inset, the time dependence of the on- 
site densities is shown. Right panel - Effective-Hamiltonian evolution of the average density within a cluster of different size where the dimensionless 
time is defined as = %~t. We show results for N = 72 and 96 sites with open-boundary conditions, obtained by the time-evolving block decimation 
technique 17 . We note that even small clusters of doubly occupied sites can effectively freeze the dynamical evolution in the large- U regime. 



an on-site energy €i = {i\H\i) and is coupled only to its nearest 
neighbors by hopping elements £*_».,■+ 1 and x 22 . It is easy to 
realize that the unitary evolution of the original many body problem 
is thus fully equivalent to the dynamics of a single particle, initially 
sitting at site 0, that is then let propagate along such a tight-binding 
chain of many-body states. We note that, both and £*_*,■+ 1 largely 
fluctuate from site to site, therefore resembling an effective Anderson 
model, even though those parameters are in reality deterministic, see 
Fig. 3. In the same figure, we also show the mean distance traveled by 
the particle after time t starting from the first site of the chain, which 
corresponds to an initial state |0) = (. . . 2, 0, 2, 0,. . .), and for different 
U/J. We observe that, for small values of (7//, the particle diffuses and 
its wave-packet finally spreads over the whole chain, in a rather 
uniform way. On the contrary, above a certain critical value of the 
interaction (U /J)f n y the particle stays localized near the origin for 
arbitrarily long times. A closer look to the structure of the on site 
energies reveal the existence of a potential well at the edge of the 
chain. This is crucial in order to understand the observed localization 
transition. Since the potential well cannot induce a true bound state 
below the bottom of the spectrum, at most a resonance will form in 
the spectrum, from which the particle could in principle escape in a 
finite time. This indeed happens at small U but apparently not at 
large (7, where the increased depth of the well and the effective 
randomness of the on-site energies conspire together to keep the 
particle localized close to the edge, preventing the the states in the 



well to hybridize with other states along the chain. We now see how 
this result connects with the previous analysis on the density relaxa- 
tion times. In the small U/J regime the particle is able to escape from 
the well and to explore larger portions of the chain, thus resulting 
into a fast density relaxation. As opposite, for large U/J, the particle 
bounces back and forth inside the well, finding hard time to escape 
from it. This lack of diffusion results into a very long-time scale for 
the density to relax to its homogeneous value. 

The above results show explicitly that some kind of localization in 
the many-body configurational space does occur, at least in the finite 
system 3 ' 23 . While such an intriguing behavior might well be a subtle 
effect due to the finite size spectrum, it could also signal the onset of a 
genuine localization that survives in the thermodynamic limit. 

We conclude the discussion on the inhomogeneous initial states 
by studying the case at density n<l, when at equilibrium there is no 
longer a Mott transition. For example, we consider n = 2/3 and 
an initial density profile (...2, 0, 0, 2, 0, 0,...). Interestingly, we find 
quite a different behavior for the density relaxation times t r (see 
Supplementary Material), with a much smoother crossover from 
small to large values of U/J and no evidence of any increase in the 
relaxation times with the system size. This fact suggests a deeper 
connection between the observed dynamical behavior and the zero 
temperature Mott transition that occurs at equilibrium and at in- 
teger filling, as suggested by calculations with infinite-coordination 
lattices 24 ' 25 . 




0 100 200 300 400 500 0 KM) 200 300 400 500 0 100 200 300 400 SOO 



site site t 

Figure 3 | Left panels - On site energies and nearest-neighbor hoppings of the effective chain that represents the Hamiltonian in the Lanczos basis 
starting from the state (...,2, 0,2,0, ...). Red points refer to U= 2/, when the particle does diffuse starting from site 0, while blue points to U= 10/, when 
it does not. The shaded regions correspond to energies less or equal that of the initial state. Right panel- Time-dependent expectation value of the wave- 
packet position of the effective particle traveling in the Hilbert space generated by a chain of R max = 1 000 Lanczos states. The red points correspond to U = 
2 / and the blue points to U = 10/ (the original lattice size of the Bose-Hubbard model is N = 12). The shaded region marks the center of the Lanczos chain, 
which is not reached in the localized regime. 
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Homogeneous initial states and Quantum Quenches. In light of 
the previous results, one may question that by choosing an 
inhomogeneous configuration of doublons we pick up a rather 
specific initial state in the Hilbert space. We are now going to 
show that the above findings strongly affect the dynamics starting 
from a perfectly homogeneous state. In this respect, a particularly 
interesting class of initial states are the ground states of H for given 
values of the interaction U b which are let evolve under the 
Hamiltonian dynamics after a sudden increase of the interaction to 
a final value Uf > U b the so-called quantum quench. Kollath and 
coworkers 26 reported evidence for the existence of two separated 
regimes in which either thermal or non-thermal behavior is 
observed for local observables. The origin of the non-thermal 
behavior in the large Uf region and the possibility of an ergodicity 
breaking in the thermodynamic limit is still highly debated 29 ' 30 . 

In the following, we focus on an average density n = 1 and we 
show that signatures of long lived metastable states of doblons can be 
identified in the dynamics after a quantum quench. At variance with 
the previous numerical calculations, now both the initial state and 
the quantum Hamiltonian do preserve the spatial homogeneity and, 
therefore, the quest for possible signatures of ergodicity breaking 
requires a different approach. Since we have identified density 
relaxation as the slowest process in the problem, we monitor the 
dynamics of the system by measuring the auto -correlation of the 
density averaged over all sites, namely through 

c(f)=^E( M 'W"'(°))-^w)^( 0 ))- ( 3 ) 

i 

For any finite size system, lim f _ >00 C(t) = 0 since the densities dec- 
orrelate at very long times. Indeed, as shown in Fig. 4, for small Uf<£i 
Uf this quantity has a very fast transient to zero. On the contrary, for 
Uf^> Uj the density auto-correlation C(t) gets stuck into a long-lived 
finite value plateau C* # 0, before approaching zero only on a much 
longer time scale. If we extract a relaxation time from C(t), we find a 
similar behavior as in Fig. 2, i.e., a dramatic increase of the relaxation 
times above a threshold value of the final interaction strength. 

In agreement with the previous analysis, the appearance of such a 
long-lived metastable state characterized by a finite plateau C* of the 
density auto -correlation function might indicate an excess of double 
occupancies that have no channel to relax. In other words, the 
dynamical constraints brought by the interaction severely slow down 
density excitations, whose characteristic time scales increase abruptly 
after a critical threshold. The main phenomenological traits of this 
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Figure 4 | Inverse relaxation times of density excitations in the 
homogeneous system, see Eq. (3). From left to right, different curves 
correspond to different initial states at UJJ= 0, 1, and 2. Insets: real part of 
the density correlations C(r) in the ergodic and in the non-ergodic region. 
Data are obtained with exact- diagonalization on a lattice with N = 12. 



dynamical arrest characterized by long-lived inhomogeneous states 
closely remind the physics of glassy materials. 

Variational description, lack of thermalization and higher dimen- 
sions. The previous discussion revealed the existence of a threshold 
energy above which a steep increase of the relaxation time of density 
fluctuations takes place. In order to assess the relevance of this 
phenomenon for the dynamics of larger systems or even for higher 
spatial dimensions, it is desirable to devise a comprehensive 
alternative framework able to catch its very characteristics. Here, 
we introduce an approach based on real-time variational Monte 
Carlo (see Supplementary Material for details), which has two 
important advantages: it allows us to follow the evolution for times 
comparable to those accessible experimentally, which are much 
longer than t-DMRG 31 ' 32 ; it can be easily extended to higher 
dimensions. We mention that mean-field-like variational ap- 
proaches to the real-time dynamics of correlated systems have 
been developed in recent years 24,25 . 

However, although they seem to capture well the main features 
of the dynamical evolution, these methods are unable to describe 
important aspects such as damping and relaxation. On the contrary, 
our approach is sufficiently rich to account for damping and relaxa- 
tion of local observables. It is based on a very simple and transparent 
out- of- equilibrium extension of the Jastrow-like variational wave 
function that was shown to describe quite accurately the equilibrium 
phase diagram of the Bose- Hubbard model 33 : 

|Y(f))=exp^Vy(n^;f)W), (4) 

where l^o) is the initial state and V^n,-, rip t) is a Jastrow factor that 
depends on the occupancies n t and rij of two sites i and; and varies 
with time so to maintain the time evolution as close as possible to the 
true evolution via the Schroedinger equation (see Supplementary 
Material). The comparison between our approach and the t- 
DMRG 34 is reported in the Supplementary Materials, demonstrating 
the high accuracy of the time-dependent variational Monte Carlo. 

Results for the time evolution of a local observable, such as the 
potential energy, after a sudden quench from U { = 2 J to a final Uf are 
shown in Fig. 5. The values of the thermal averages have been com- 
puted in the grand- canonical ensemble by means of finite - 
temperature quantum Monte Carlo calculations 35 , with the effective 
temperature fixed by the average energy of the initial state, which we 
take as the best variational approximation for the ground state at U { 
= 2J. 

As shown in Fig. 5, in the region of small Uf we observe a damping 
of the average potential- energy, which approaches a quasi- steady 
stationary value in contrast to the simple Gutzwiller wave func- 
tion 24 ' 25 . In this regime, damping is mainly due to a density- density 
Jastrow factor of the form Vifn b rip) = Vy(f) n t rip which already at 
equilibrium was shown to provide a satisfactory description of the 
physical behavior 33 . This fact enlightens the relevance of the 
Bogoliubov modes whose dephasing during the time evolution 
allows to approach the stationary state. Remarkably, the steady state 
averages coincide with the thermal ones; a signal that the dynamics is 
ergodic. 

In the region of large interactions Up a simple density- density 
Jastrow factor does not account for all relaxation pathways, which 
will now mainly result from specific correlations among doublons, 
holons and between holons and doublons. The effective Hamiltonian 
(2) indeed explicitly shows that doublons attract each other as well as 
holons do, while doublons repel holons. These correlations, as well as 
other among higher on-site densities, can be easily implemented via 
the the Jastrow factor in Eq. (4) and indeed substantially improve the 
dynamics. Interestingly, the effective interaction between doblons 
that results from the dynamical variational calculation turns to be 
attractive, therefore leading to a consistent determination of the 
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Figure 5 | One-dimensional results for the time- dependent expectation 
values of the on-site potential energy s p ( t) = ^ (n,- ( n t — 1) ) in the ergodic 
(upper panel) and in the non-ergodic regions (lower panel). The initial 
state is the ground state of the Bose-Hubbard Hamiltonian with U r ■ = 2 J 
and the considered system size is N = 200. Grand-canonical thermal 
averages are shown for comparison as dashed horizontal lines. 



high- energy incoherent excitations that do not have channels to relax 
efficiently. This scenario is also consistent with the observed strong 
doping dependence of the relaxation time (see Supplementary 
Material). 

Discussion 

In conclusion, we have found that the dynamical constraints brought 
by a strong interaction can trap the evolution of repulsive bosons 
hopping on a lattice into metastable states that lack translational 
symmetry, provided that the energy stored into the initial state is 
above a threshold. We pointed out that a self-induced effective 
attraction among doublons is one of the major processes that can 
effectively freeze the dynamics on long time scales. Such a mech- 
anism is recognized to play a role in the density relaxation processes 
of purely homogeneous systems through a dynamical arrest visible in 
time-dependent density correlations. The main features of this intri- 
guing behavior, namely the slowing down of density excitations and 
the long-lived inhomogeneous pattern, resembles closely a kind of 
glass transition. 

Moreover, we have shown that the time evolution of the many- 
body problem can be mapped onto that of a particle moving from the 
edge of a semi-infinite tight-binding chain with nearest-neighbor 
hopping, where each site represents a many-body wave function. 
This model looks like an Anderson model, since both the on-site 
energy and the hopping vary from site to site, with a potential well 
at one edge due to the high-energy content of the initial state. 
Interestingly, we find a delocalization-localization transition in this 
problem, with the particle being unable to diffuse on the whole chain 
above a certain value of the well depth. We consider this analogy 
quite suggestive and potentially constitute an even stronger indica- 
tion of ergodicity breaking in the many-body space which is worth to 
be further investigated. 



anticipated dynamical effects that drive the dynamics in this regime. 
As we see from Fig. 5, in the region of very large Uf the potential- 
energy expectation values do show a damping to a non -thermal 
quasi-steady state on a time of the order of t d ~ 1/Uf. This fast time 
scale must be put in comparison with the much longer one, t r of 
Fig. 4. It is natural to identify t r with the time scale that controls the 
eventual escape from the quasi- steady state, hence the approach to 
thermal equilibrium. Whether this time scale does truly diverge in 
the thermodynamic limit, or rather saturates to a very large value 
which is still out of reach for state of the art numerics, it is certainly an 
important issue that cannot be definitively solved. However, we can 
safely state that a large finite system of actual experimental relevance 
will get stuck for a quite long time into highly inhomogeneous meta- 
stable states, which we revealed to be on the verge of a spatial sym- 
metry breaking. 

The above results have been obtained for out-of-equilibrium one 
dimensional systems and therefore leave the questions concerning 
the dependence on the dimensionality of the problem still open. 
However, the anomalously long-time relaxation of the density 
auto- correlation points towards a kind of glassy behavior that should 
be observable even in higher dimensions, provided that the inter- 
action induces sufficiently strong dynamical constraints. In this 
regard, we have studied the two-dimensional case by means of our 
time-dependent variational scheme, verifying that a similar behavior 
occurs even in two dimensions. In Fig. 6, we show the results of the 
potential-energy expectation values as a function of the time. As 
before, if the repulsion is weak, the time average coincides with the 
thermal ones, while for strong repulsion there is a clear difference 
between thermal and time averages, giving rise to a nonergodic 
dynamics. We therefore argue that the phenomenology we have 
hereby identified is almost independent on the dimensionality and 
it is rather due to the existence in strongly correlated systems of 
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Figure 6 | Two dimensional results for the time-dependent expectation 
values of the on-site potential energy s p ( t) = (tii (ti{ — 1 ) ) in the ergodic 
(upper panel) and in the non-ergodic regions (lower panel). The initial 
state is the ground state of the Bose-Hubbard Hamiltonian with U t ■= 4 J 
and the considered system size is N = 20X20. Grand-canonical thermal 
averages are shown for comparison as dashed horizontal lines. 
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